
## Source Code

source("NewSim.R")
library(survival)
set.seed(3)
rxe<-(-1)
base<-1
maxev<-10
maxt <- 100
p<-1

# Vary censoring rates from 0 to 50% by 10% steps
rates <- (0:5)/10

## Monte Carlo Analysis

nn <- 100 #Level-2 sample size
nrep <-1000 #number of repetitions

####################################
##      EXPONENTIAL
####################################


# Event Dependence and Heterogeneity

for (rate in rates){
  eventcount <- c(rep(NA,nrep))
  censrate <- c(rep(NA,nrep))
  samplesize <- c(rep(NA,nrep))
  
  betaag <- c(rep(NA,nrep))
  stderag <-  c(rep(NA,nrep))
  naiveag <-  c(rep(NA,nrep))

  betacondelapse <-  c(rep(NA,nrep))
  stdercondelapse <-  c(rep(NA,nrep))
  naivecondelapse <-  c(rep(NA,nrep))

  betacondgap <- c(rep(NA,nrep))
  stdercondgap <- c(rep(NA,nrep))
  naivecondgap <- c(rep(NA,nrep))

  betafrailelapse <- c(rep(NA,nrep))
  stderfrailelapse <- c(rep(NA,nrep))
  thetafrailelapse <- c(rep(NA,nrep))
  thetapvalfrailelapse <- c(rep(NA,nrep))

  betafrailgap <-  c(rep(NA,nrep))
  stderfrailgap <-  c(rep(NA,nrep))
  thetafrailgap <- c(rep(NA,nrep))
  thetapvalfrailgap <- c(rep(NA,nrep))

  cat("Censoring Rate:\t", rate, "\n")
  for (i in 1:nrep){
    set.seed(3*i)
    if (as.integer(i/100)==(i/100)) { cat("Iterations = ", i, "\n") }
    sim1<-simdatashapegamma.random(nn,1,maxev,rate, maxt, p, rxe, base, klambdaweibull, 0, 1)
    
                                        # Cond, Frailty Gap
    frailgap<-coxph(formula = Surv(sim1[,9],sim1[,10],sim1[,4],type = "counting")~ sim1[,5] + strata(sim1[,7]) + frailty(sim1[,1]), na.action = na.exclude, eps = 0.0001, iter.max = 100, method = "efron")
    betafrailgap[i]<-frailgap$coef
    stderfrailgap[i]<-sqrt(diag(frailgap$var))
    thetafrailgap[i] <- frailgap$history[[1]]$history[length(frailgap$history[[1]]$history[,1]),1]
    thetapvalfrailgap[i] <- (frailgap$printfun[[1]])(frailgap$frail, frailgap$fvar,, frailgap$df[2], frailgap$history[[1]])$coef[6]
    
                                        # Frailty, Elasped
    frailelapse<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4],type = "counting")~ sim1[,5] + frailty(sim1[,1]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron")
    betafrailelapse[i] <- frailelapse$coef
    if(is.na(frailelapse$coef)){ stderfrailelapse[i]<-NA}  else { stderfrailelapse[i]<-sqrt(diag(frailelapse$var)) }
    if(is.na(frailelapse$coef)){ thetafrailelapse[i]<-NA} else {thetafrailelapse[i] <-  frailelapse$history[[1]]$history[length(frailelapse$history[[1]]$history[,1]),1]}
    if(is.na(frailelapse$coef)){ thetapvalfrailelapse[i]<-NA} else {thetapvalfrailelapse[i] <- (frailelapse$printfun[[1]])(frailelapse$frail, frailelapse$fvar,, frailelapse$df[2], frailelapse$history[[1]])$coef[6]}
    
  # Andersen-Gill
    agevent<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4], type = "counting")~ sim1[,5] + cluster(sim1[,1]), na.action = na.exclude, eps =  0.0001, iter.max = 10, method = "efron", robust = T)
    betaag[i]<-agevent$coef
    stderag[i]<-sqrt(diag(agevent$var))
    naiveag[i]<-sqrt(diag(agevent$naive.var))
    
                                        # Conditional Gap
    condgap<-coxph(formula = Surv(sim1[,9],sim1[,10],sim1[,4],type = "counting")~ sim1[,5] + cluster(sim1[,1]) + strata(sim1[,7]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron", robust = T)
  betacondgap[i] <- condgap$coef
    stdercondgap[i]<-sqrt(diag(condgap$var))
    naivecondgap[i]<-sqrt(diag(condgap$naive.var))
  
  # Conditional Elapsed
    condelapse<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4],type = "counting")~ sim1[,5] + cluster(sim1[,1]) + strata(sim1[,7]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron", robust = T)
    betacondelapse[i] <-condelapse$coef
    stdercondelapse[i]<-sqrt(diag(condelapse$var))
    naivecondelapse[i]<-sqrt(diag(condelapse$naive.var))

  # Event Count
    zen <- aggregate(x=sim1$enum, by=list(sim1$id), FUN="max")
    eventcount[i] <- mean(zen$x)
    censrate[i] <- 1-mean(sim1$status)
    samplesize[i] <- nrow(sim1)
  }
  simnumber <- c(1:nrep)
  simdata <- cbind(simnumber, betafrailgap, stderfrailgap, thetafrailgap, thetapvalfrailgap, betafrailelapse, stderfrailelapse, thetafrailelapse, thetapvalfrailelapse, betaag, stderag, naiveag, betacondgap, stdercondgap, naivecondgap, betacondelapse, stdercondelapse, naivecondelapse, eventcount, censrate, samplesize)
  filestring <- paste("exp_ed_het_randomcens", rate, ".txt", sep="")
  write.table(simdata, file = filestring, sep="\t", row.names=FALSE)
}


## No Event Dependence, No Heterogeneity


nn <- 100


for (rate in rates){
  eventcount <- c(rep(NA,nrep))
  censrate <- c(rep(NA,nrep))
  samplesize <- c(rep(NA,nrep))
  
  betaag <- c(rep(NA,nrep))
  stderag <-  c(rep(NA,nrep))
  naiveag <-  c(rep(NA,nrep))
  
  betacondelapse <-  c(rep(NA,nrep))
  stdercondelapse <-  c(rep(NA,nrep))
  naivecondelapse <-  c(rep(NA,nrep))

  betacondgap <- c(rep(NA,nrep))
  stdercondgap <- c(rep(NA,nrep))
  naivecondgap <- c(rep(NA,nrep))
  
  betafrailelapse <- c(rep(NA,nrep))
  stderfrailelapse <- c(rep(NA,nrep))
  thetafrailelapse <- c(rep(NA,nrep))
  thetapvalfrailelapse <- c(rep(NA,nrep))

  betafrailgap <-  c(rep(NA,nrep))
  stderfrailgap <-  c(rep(NA,nrep))
  thetafrailgap <- c(rep(NA,nrep))
  thetapvalfrailgap <- c(rep(NA,nrep))

  cat("Censoring Rate:\t", rate, "\n")
  for (i in 1:nrep){
    
    if (as.integer(i/100)==(i/100)) { cat("Iterations = ", i, "\n") }
    set.seed(3*i)
    sim1<-simdatashapegamma.random(nn,1,maxev,rate, maxt, p, rxe, base, klambdaweibull, 1, 0)
    
  # Cond, Frailty Gap
    frailgap<-coxph(formula = Surv(sim1[,9],sim1[,10],sim1[,4],type = "counting")~ sim1[,5] + strata(sim1[,7]) + frailty(sim1[,1]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron")
    betafrailgap[i]<-frailgap$coef
    stderfrailgap[i]<-sqrt(diag(frailgap$var))
    thetafrailgap[i] <- frailgap$history[[1]]$history[length(frailgap$history[[1]]$history[,1]),1]
    thetapvalfrailgap[i] <- (frailgap$printfun[[1]])(frailgap$frail, frailgap$fvar,, frailgap$df[2], frailgap$history[[1]])$coef[6]
    
  # Frailty, Elasped
    frailelapse<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4],type = "counting")~ sim1[,5] + frailty(sim1[,1]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron")
    betafrailelapse[i] <- frailelapse$coef
    stderfrailelapse[i]<-sqrt(diag(frailelapse$var))
    thetafrailelapse[i] <-  frailelapse$history[[1]]$history[length(frailelapse$history[[1]]$history[,1]),1]
    thetapvalfrailelapse[i] <- (frailelapse$printfun[[1]])(frailelapse$frail, frailelapse$fvar,, frailelapse$df[2], frailelapse$history[[1]])$coef[6]
    
  # Andersen-Gill
    agevent<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4], type = "counting")~ sim1[,5] + cluster(sim1[,1]), na.action = na.exclude, eps =  0.0001, iter.max = 10, method = "efron", robust = T)
    betaag[i]<-agevent$coef
    stderag[i]<-sqrt(diag(agevent$var))
    naiveag[i]<-sqrt(diag(agevent$naive.var))
  
  # Conditional Gap
    condgap<-coxph(formula = Surv(sim1[,9],sim1[,10],sim1[,4],type = "counting")~ sim1[,5] + cluster(sim1[,1]) + strata(sim1[,7]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron", robust = T)
    betacondgap[i] <- condgap$coef
    stdercondgap[i]<-sqrt(diag(condgap$var))
    naivecondgap[i]<-sqrt(diag(condgap$naive.var))
  
  # Conditional Elapsed
    condelapse<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4],type = "counting")~ sim1[,5] + cluster(sim1[,1]) + strata(sim1[,7]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron", robust = T)
    betacondelapse[i] <-condelapse$coef
    stdercondelapse[i]<-sqrt(diag(condelapse$var))
    naivecondelapse[i]<-sqrt(diag(condelapse$naive.var))

  # Data summaries, Event Count per Id, censor rate, overall n
    zen <- aggregate(x=sim1$enum, by=list(sim1$id), FUN="max")
    eventcount[i] <- mean(zen$x)
    censrate[i] <- 1-mean(sim1$status)
    samplesize[i] <- nrow(sim1)
  }
  
  simnumber <- c(1:nrep)

  simdata <- cbind(simnumber, betafrailgap, stderfrailgap, thetafrailgap, thetapvalfrailgap, betafrailelapse, stderfrailelapse, thetafrailelapse, thetapvalfrailelapse, betaag, stderag, naiveag, betacondgap, stdercondgap, naivecondgap, betacondelapse, stdercondelapse, naivecondelapse, eventcount, censrate, samplesize)
  filestring <- paste("exp_noed_nohet_randomcens", rate, ".txt", sep="")
  write.table(simdata, file = filestring, sep="\t", row.names=FALSE)
}

# Event Dependence, No Heterogeneity

nn <- 100

for (rate in rates){
  eventcount <- c(rep(NA,nrep))
  censrate <- c(rep(NA,nrep))
  samplesize <- c(rep(NA,nrep))
  
  betaag <- c(rep(NA,nrep))
  stderag <-  c(rep(NA,nrep))
  naiveag <-  c(rep(NA,nrep))

  betacondelapse <-  c(rep(NA,nrep))
  stdercondelapse <-  c(rep(NA,nrep))
  naivecondelapse <-  c(rep(NA,nrep))

  betacondgap <- c(rep(NA,nrep))
  stdercondgap <- c(rep(NA,nrep))
  naivecondgap <- c(rep(NA,nrep))

  betafrailelapse <- c(rep(NA,nrep))
  stderfrailelapse <- c(rep(NA,nrep))
  thetafrailelapse <- c(rep(NA,nrep))
  thetapvalfrailelapse <- c(rep(NA,nrep))

  betafrailgap <-  c(rep(NA,nrep))
  stderfrailgap <-  c(rep(NA,nrep))
  thetafrailgap <- c(rep(NA,nrep))
  thetapvalfrailgap <- c(rep(NA,nrep))

  cat("Censoring Rate:\t", rate, "\n")
  for (i in 1:nrep){

    if (as.integer(i/100)==(i/100)) { cat("Iterations = ", i, "\n") }
    set.seed(3*i)
    sim1<-simdatashapegamma.random(nn,1,maxev,rate, maxt, p, rxe, base, klambdaweibull, 0, 0)
    
  # Cond, Frailty Gap
    frailgap<-coxph(formula = Surv(sim1[,9],sim1[,10],sim1[,4],type = "counting")~ sim1[,5] + strata(sim1[,7]) + frailty(sim1[,1]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron")
    betafrailgap[i]<-frailgap$coef
    stderfrailgap[i]<-sqrt(diag(frailgap$var))
    thetafrailgap[i] <- frailgap$history[[1]]$history[length(frailgap$history[[1]]$history[,1]),1]
    thetapvalfrailgap[i] <- (frailgap$printfun[[1]])(frailgap$frail, frailgap$fvar,, frailgap$df[2], frailgap$history[[1]])$coef[6]
    
  # Frailty, Elasped
    frailelapse<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4],type = "counting")~ sim1[,5] + frailty(sim1[,1]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron")
    betafrailelapse[i] <- frailelapse$coef
    stderfrailelapse[i]<-sqrt(diag(frailelapse$var))
    thetafrailelapse[i] <-  frailelapse$history[[1]]$history[length(frailelapse$history[[1]]$history[,1]),1]
    thetapvalfrailelapse[i] <- (frailelapse$printfun[[1]])(frailelapse$frail, frailelapse$fvar,, frailelapse$df[2], frailelapse$history[[1]])$coef[6]
    
  # Andersen-Gill
    agevent<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4], type = "counting")~ sim1[,5] + cluster(sim1[,1]), na.action = na.exclude, eps =  0.0001, iter.max = 10, method = "efron", robust = T)
    betaag[i]<-agevent$coef
    stderag[i]<-sqrt(diag(agevent$var))
    naiveag[i]<-sqrt(diag(agevent$naive.var))
  
  # Conditional Gap
    condgap<-coxph(formula = Surv(sim1[,9],sim1[,10],sim1[,4],type = "counting")~ sim1[,5] + cluster(sim1[,1]) + strata(sim1[,7]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron", robust = T)
    betacondgap[i] <- condgap$coef
    stdercondgap[i]<-sqrt(diag(condgap$var))
    naivecondgap[i]<-sqrt(diag(condgap$naive.var))
  
  # Conditional Elapsed
    condelapse<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4],type = "counting")~ sim1[,5] + cluster(sim1[,1]) + strata(sim1[,7]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron", robust = T)
    betacondelapse[i] <-condelapse$coef
    stdercondelapse[i]<-sqrt(diag(condelapse$var))
    naivecondelapse[i]<-sqrt(diag(condelapse$naive.var))

  # Data summaries, Event Count per Id, censor rate, overall n
    zen <- aggregate(x=sim1$enum, by=list(sim1$id), FUN="max")
    eventcount[i] <- mean(zen$x)
    censrate[i] <- 1-mean(sim1$status)
    samplesize[i] <- nrow(sim1)
  }
  simnumber <- c(1:nrep)

  simdata <- cbind(simnumber, betafrailgap, stderfrailgap, thetafrailgap, thetapvalfrailgap, betafrailelapse, stderfrailelapse, thetafrailelapse, thetapvalfrailelapse, betaag, stderag, naiveag, betacondgap, stdercondgap, naivecondgap, betacondelapse, stdercondelapse, naivecondelapse, eventcount, censrate, samplesize)
  filestring <- paste("exp_ed_nohet_randomcens", rate, ".txt", sep="")
  write.table(simdata, file = filestring, sep="\t", row.names=FALSE)
}


# No Event Dependence, Heterogeneity

nn<-100

for (rate in rates){
  eventcount <- c(rep(NA,nrep))
  censrate <- c(rep(NA,nrep))
  samplesize <- c(rep(NA,nrep))
  
  betaag <- c(rep(NA,nrep))
  stderag <-  c(rep(NA,nrep))
  naiveag <-  c(rep(NA,nrep))

  betacondelapse <-  c(rep(NA,nrep))
  stdercondelapse <-  c(rep(NA,nrep))
  naivecondelapse <-  c(rep(NA,nrep))

  betacondgap <- c(rep(NA,nrep))
  stdercondgap <- c(rep(NA,nrep))
  naivecondgap <- c(rep(NA,nrep))

  betafrailelapse <- c(rep(NA,nrep))
  stderfrailelapse <- c(rep(NA,nrep))
  thetafrailelapse <- c(rep(NA,nrep))
  thetapvalfrailelapse <- c(rep(NA,nrep))
  
  betafrailgap <-  c(rep(NA,nrep))
  stderfrailgap <-  c(rep(NA,nrep))
  thetafrailgap <- c(rep(NA,nrep))
  thetapvalfrailgap <- c(rep(NA,nrep))
  

  cat("Censoring Rate:\t", rate, "\n")
  for (i in 1:nrep){
    
    if (as.integer(i/100)==(i/100)) { cat("Iterations = ", i, "\n") }
    set.seed(3*i)
    sim1<-simdatashapegamma.random(nn,1,maxev,rate, maxt, p, rxe, base, klambdaweibull, 1, 1)
    
    # Cond, Frailty Gap
    frailgap<-coxph(formula = Surv(sim1[,9],sim1[,10],sim1[,4],type = "counting")~ sim1[,5] + strata(sim1[,7]) + frailty(sim1[,1]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron")
    betafrailgap[i]<-frailgap$coef
    stderfrailgap[i]<-sqrt(diag(frailgap$var))
    thetafrailgap[i] <- frailgap$history[[1]]$history[length(frailgap$history[[1]]$history[,1]),1]
    thetapvalfrailgap[i] <- (frailgap$printfun[[1]])(frailgap$frail, frailgap$fvar,, frailgap$df[2], frailgap$history[[1]])$coef[6]
    
  # Frailty, Elasped
    frailelapse<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4],type = "counting")~ sim1[,5] + frailty(sim1[,1]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron")
    betafrailelapse[i] <- frailelapse$coef
    stderfrailelapse[i]<-sqrt(diag(frailelapse$var))
    thetafrailelapse[i] <-  frailelapse$history[[1]]$history[length(frailelapse$history[[1]]$history[,1]),1]
    thetapvalfrailelapse[i] <- (frailelapse$printfun[[1]])(frailelapse$frail, frailelapse$fvar,, frailelapse$df[2], frailelapse$history[[1]])$coef[6]
    
  # Andersen-Gill
    agevent<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4], type = "counting")~ sim1[,5] + cluster(sim1[,1]), na.action = na.exclude, eps =  0.0001, iter.max = 10, method = "efron", robust = T)
    betaag[i]<-agevent$coef
    stderag[i]<-sqrt(diag(agevent$var))
    naiveag[i]<-sqrt(diag(agevent$naive.var))
    
  # Conditional Gap
    condgap<-coxph(formula = Surv(sim1[,9],sim1[,10],sim1[,4],type = "counting")~ sim1[,5] + cluster(sim1[,1]) + strata(sim1[,7]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron", robust = T)
    betacondgap[i] <- condgap$coef
    stdercondgap[i]<-sqrt(diag(condgap$var))
    naivecondgap[i]<-sqrt(diag(condgap$naive.var))
  
  # Conditional Elapsed
    condelapse<-coxph(formula = Surv(sim1[,2],sim1[,3],sim1[,4],type = "counting")~ sim1[,5] + cluster(sim1[,1]) + strata(sim1[,7]), na.action = na.exclude, eps = 0.0001, iter.max = 10, method = "efron", robust = T)
    betacondelapse[i] <-condelapse$coef
    stdercondelapse[i]<-sqrt(diag(condelapse$var))
    naivecondelapse[i]<-sqrt(diag(condelapse$naive.var))

  # Event Count
    zen <- aggregate(x=sim1$enum, by=list(sim1$id), FUN="max")
    eventcount[i] <- mean(zen$x)
    censrate[i] <- 1-mean(sim1$status)
    samplesize[i] <- nrow(sim1)
  }
  
  simnumber <- c(1:nrep)
  
  simdata <- cbind(simnumber, betafrailgap, stderfrailgap, thetafrailgap, thetapvalfrailgap, betafrailelapse, stderfrailelapse, thetafrailelapse, thetapvalfrailelapse, betaag, stderag, naiveag, betacondgap, stdercondgap, naivecondgap, betacondelapse, stdercondelapse, naivecondelapse, eventcount, censrate, samplesize)
  filestring <- paste("exp_noed_het_randomcens", rate, ".txt", sep="")
  write.table(simdata, file = filestring, sep="\t", row.names=FALSE)
}

